% find playback onset cells
low_t=5;
high_t=20;
snip_length=low_t+high_t+1;
res=[];
cnt2=1;
for lnd=1:4
    for fnd=1:size(proj_meta(cur_site).raw_data(lnd,1).smoothed_ROIs,2)
        snips={};
        for tnd=1:length(proj_meta(cur_site).ExpGroup)
            raw=proj_meta(cur_site).raw_data(lnd,tnd).smoothed_ROIs(fnd).activity;
            raw_act=raw/median(raw);
            cnt=1;
            ons=find(diff(proj_meta(cur_site).raw_data(lnd,tnd).velP_ind)==1);
            for gnd=1:length(ons)
                if ons(gnd)-low_t>0 && ons(gnd)+high_t<length(proj_meta(cur_site).raw_data(lnd,tnd).velP_ind)
                    cur_low=ons(gnd)-low_t:ons(gnd);
                    cur_high=ons(gnd)+1:ons(gnd)+high_t;
                    if sum(proj_meta(cur_site).raw_data(lnd,tnd).velP_ind(cur_low))==0 && sum(proj_meta(cur_site).raw_data(lnd,tnd).velP_ind(cur_high))==high_t
                        if sum(proj_meta(cur_site).raw_data(lnd,tnd).velM_ind(cur_low:cur_high))==0
                            snips{tnd}(cnt,:)=raw_act([cur_low cur_high]);
                            cnt=cnt+1;
                        end
                    end
                end
            end
        end
        exp_av=[];
        for tnd=1:length(proj_meta(cur_site).ExpGroup)
            if ~isempty(snips{tnd})
                exp_av(tnd,:)=mean(snips{tnd},1);
            else
                exp_av(tnd,:)=0;
            end
        end
        if sum(mean(exp_av(:,1:3),2)*1.1<mean(exp_av(:,snip_length-10:snip_length-8),2))
            res(cnt2,1)=lnd;
            res(cnt2,2)=fnd;
            cnt2=cnt2+1;
        end
    end
end
vis_res=res;
all_res={};
for nnd=1:size(vis_res,1)
    snips={};
    for tnd=1:length(proj_meta(cur_site).ExpGroup)
        raw=proj_meta(cur_site).raw_data(vis_res(nnd,1),tnd).smoothed_ROIs(vis_res(nnd,2)).activity;
        raw_act=raw/median(raw);
        cnt=1;
        ons=find(diff(proj_meta(cur_site).raw_data(vis_res(nnd,1),tnd).velP_ind)==1);
        for gnd=1:length(ons)
            if ons(gnd)-low_t>0 && ons(gnd)+high_t<length(proj_meta(cur_site).raw_data(vis_res(nnd,1),tnd).velP_ind)
                cur_low=ons(gnd)-low_t:ons(gnd);
                cur_high=ons(gnd)+1:ons(gnd)+high_t;
                if sum(proj_meta(cur_site).raw_data(vis_res(nnd,1),tnd).velP_ind(cur_low))==0 && sum(proj_meta(cur_site).raw_data(vis_res(nnd,1),tnd).velP_ind(cur_high))==high_t
                    if sum(proj_meta(cur_site).raw_data(vis_res(nnd,1),tnd).velM_ind(cur_low:cur_high))==0
                        snips{tnd}(cnt,:)=raw_act([cur_low cur_high]);
                        cnt=cnt+1;
                    end
                end
            end
        end
    end
    for tnd=1:length(proj_meta(cur_site).ExpGroup)
        all_res{tnd}(nnd,:)=mean(snips{tnd},1);
    end
end
vis_all_res=all_res;
figure;
for tnd=1:length(proj_meta(cur_site).ExpGroup)
    set(gcf,'position',[2003         649        1800         336])
    subplot(1,length(proj_meta(cur_site).ExpGroup),tnd)
    imagesc(vis_all_res{tnd})
    set(gca,'clim',[0.95 1.85])
end
title('Playback onset')
drawnow

%find running onset cells
low_t=5;
high_t=20;
snip_length=low_t+high_t+1;
res=[];
cnt2=1;
for lnd=1:4
    for fnd=1:size(proj_meta(cur_site).raw_data(lnd,1).smoothed_ROIs,2)
        snips={};
        for tnd=1:length(proj_meta(cur_site).ExpGroup)
            raw=proj_meta(cur_site).raw_data(lnd,tnd).smoothed_ROIs(fnd).activity;
            raw_act=raw/median(raw);
            cnt=1;
            ons=find(diff(proj_meta(cur_site).raw_data(lnd,tnd).velM_ind)==1);
            for gnd=1:length(ons)
                if ons(gnd)-low_t>0 && ons(gnd)+high_t<length(proj_meta(cur_site).raw_data(lnd,tnd).velM_ind)
                    cur_low=ons(gnd)-low_t:ons(gnd);
                    cur_high=ons(gnd)+1:ons(gnd)+high_t;
                    if sum(proj_meta(cur_site).raw_data(lnd,tnd).velM_ind(cur_low))==0 && sum(proj_meta(cur_site).raw_data(lnd,tnd).velM_ind(cur_high))==high_t
                        snips{tnd}(cnt,:)=raw_act([cur_low cur_high]);
                        cnt=cnt+1;
                    end
                end
            end
        end
        exp_av=[];
        for tnd=1:length(proj_meta(cur_site).ExpGroup)
            if ~isempty(snips{tnd})
                exp_av(tnd,:)=mean(snips{tnd},1);
            else
                exp_av(tnd,:)=0;
            end
        end
        if sum(mean(exp_av(:,1:3),2)*1.1<mean(exp_av(:,snip_length-10:snip_length-8),2))
            res(cnt2,1)=lnd;
            res(cnt2,2)=fnd;
            cnt2=cnt2+1;
        end
    end
end
run_res=res;
remove=[];
for snd=1:size(vis_res,1)
    if sum(sum(ismember(run_res,vis_res(snd,:)),2)==2)
        remove(snd)=find(sum(ismember(run_res,vis_res(snd,:)),2)==2);
    end
end
remove(remove==0)=[];
run_res(remove,:)=[];

all_res={};
for nnd=1:size(run_res,1)
    snips={};
    for tnd=1:length(proj_meta(cur_site).ExpGroup)
        raw=proj_meta(cur_site).raw_data(run_res(nnd,1),tnd).smoothed_ROIs(run_res(nnd,2)).activity;
        raw_act=raw/median(raw);
        cnt=1;
        ons=find(diff(proj_meta(cur_site).raw_data(run_res(nnd,1),tnd).velM_ind)==1);
        for gnd=1:length(ons)
            if ons(gnd)-low_t>0 && ons(gnd)+high_t<length(proj_meta(cur_site).raw_data(run_res(nnd,1),tnd).velM_ind)
                cur_low=ons(gnd)-low_t:ons(gnd);
                cur_high=ons(gnd)+1:ons(gnd)+high_t;
                if sum(proj_meta(cur_site).raw_data(run_res(nnd,1),tnd).velM_ind(cur_low))==0 && sum(proj_meta(cur_site).raw_data(run_res(nnd,1),tnd).velM_ind(cur_high))==high_t
                    snips{tnd}(cnt,:)=raw_act([cur_low cur_high]);
                    cnt=cnt+1;
                end
            end
        end
    end
    for tnd=1:length(proj_meta(cur_site).ExpGroup)
        all_res{tnd}(nnd,:)=mean(snips{tnd},1);
    end
end
run_all_res=all_res;
figure;
for tnd=1:7
    set(gcf,'position',[2003         649        1800         336])
    subplot(1,7,tnd)
    imagesc(run_all_res{tnd})
    set(gca,'clim',[0.95 1.85])
end
title('Running onset')

%find pert onset response
low_t=5;
high_t=10;
snip_length=low_t+high_t+1;
pert_all_res={};
ids=[-1:-1:-6 1:7];
for xnd=1:length(ids)
    cur_pert=ids(xnd);
    res=[];
    cnt2=1;
    for lnd=1:4
        cur_site=34+lnd;
        for fnd=1:size(proj_meta(cur_site).raw_data(1).sig_act,2)
            %         cur_cell=fnd;
            snips={};
            for tnd=1:7
                raw_act=piezo_act(lnd).activity{tnd}(fnd,:);
                %             raw_act=all_act{tnd}(fnd,:);
                cnt=1;
                %             ons=find(diff(abs(proj_meta(cur_site).raw_data(tnd).stim_id))>0);
                cur_stim_id=proj_meta(cur_site).raw_data(tnd).stim_id;
                cur_stim_id(cur_stim_id~=cur_pert)=0;
                ons=find(diff(proj_meta(cur_site).raw_data(tnd).stim_id)==cur_pert);
                for gnd=1:length(ons)
                    if ons(gnd)-low_t>0 && ons(gnd)+high_t<length(proj_meta(cur_site).raw_data(tnd).runM)
                        cur_low=ons(gnd)-low_t:ons(gnd);
                        cur_high=ons(gnd)+1:ons(gnd)+high_t;
                        if sum(proj_meta(cur_site).raw_data(tnd).runM(cur_low))==low_t+1 && sum(proj_meta(cur_site).raw_data(tnd).runM(cur_high))==high_t
                            snips{tnd}(cnt,:)=raw_act([cur_low cur_high]);
                            cnt=cnt+1;
                        end
                    end
                end
            end
            exp_av=[];
            for tnd=1:7
                if ~isempty(snips{tnd})
                    exp_av(tnd,:)=mean(snips{tnd},1);
                else
                    exp_av(tnd,:)=0;
                end
            end
            if sum(mean(exp_av(:,5:7),2)*1.1<mean(exp_av(:,snip_length-7:snip_length-5),2))
                res(cnt2,1)=lnd;
                res(cnt2,2)=fnd;
                cnt2=cnt2+1;
            end
        end
    end
    pert_res=res;
    
    remove=[];
    cnt3=1;
    for snd=1:size(vis_res,1)
        if sum(sum(ismember(pert_res,vis_res(snd,:)),2)==2)
            remove(cnt3)=find(sum(ismember(pert_res,vis_res(snd,:)),2)==2);
            cnt3=cnt3+1;
        end
    end
    for snd=1:size(run_res,1)
        if sum(sum(ismember(pert_res,run_res(snd,:)),2)==2)
            remove(cnt3)=find(sum(ismember(pert_res,run_res(snd,:)),2)==2);
            cnt3=cnt3+1;
        end
    end
    pert_res(remove,:)=[];
    
    
    all_res={};
    for nnd=1:size(pert_res,1)
        cur_site=34+pert_res(nnd,1);
        snips={};
        for tnd=1:7
            raw_act=piezo_act(pert_res(nnd,1)).activity{tnd}(pert_res(nnd,2),:);
            %             raw_act=all_act{tnd}(fnd,:);
            cnt=1;
            cur_stim_id=proj_meta(cur_site).raw_data(tnd).stim_id;
                cur_stim_id(cur_stim_id~=cur_pert)=0;
                ons=find(diff(proj_meta(cur_site).raw_data(tnd).stim_id)==cur_pert);
            for gnd=1:length(ons)
                if ons(gnd)-low_t>0 && ons(gnd)+high_t<length(proj_meta(cur_site).raw_data(tnd).runM)
                    cur_low=ons(gnd)-low_t:ons(gnd);
                    cur_high=ons(gnd)+1:ons(gnd)+high_t;
                    if sum(proj_meta(cur_site).raw_data(tnd).runM(cur_low))==low_t+1 && sum(proj_meta(cur_site).raw_data(tnd).runM(cur_high))==high_t
                        snips{tnd}(cnt,:)=raw_act([cur_low cur_high]);
                        cnt=cnt+1;
                    end
                end
            end
        end
        for tnd=1:7
            all_res{tnd}(nnd,:)=mean(snips{tnd},1);
        end
    end
    pert_all_res{xnd}=all_res;
end

for dnd=1:length(ids)
    figure;
    for tnd=1:7
        set(gcf,'position',[2003         649        1800         336])
        subplot(1,7,tnd)
        imagesc(pert_all_res{dnd}{tnd})
        set(gca,'clim',[0.8 2])
    end
    title(['pert ' num2str(ids(dnd)) ' onset'])
end


av_arr={};
for dnd=1:length(ids)
    for ynd=1:7
        if size(pert_all_res{dnd}{1},1)~=1
            av_arr{dnd}(:,ynd)=mean(pert_all_res{dnd}{ynd});
        else
            av_arr{dnd}(:,ynd)=pert_all_res{dnd}{ynd};
        end
    end
end
names={};
names{1}='t1';
names{2}='t2';
names{3}='t3';
names{4}='t4';
names{5}='t5';
names{6}='t6';
names{7}='t7';
for dnd=1:length(ids)
    figure;plot(av_arr{dnd})
    legend(names,'location','northwest');
end



































































% 
% 
% 
% 
% 
% low_t=5;
% high_t=20;
% snip_length=low_t+high_t+1;
% res=[];
% cnt2=1;
% for lnd=1:4
%     cur_site=34+lnd;
%     
%     order={};
%     exp_inds={};
%     for tnd=1:7
%         nbr_frames=proj_meta(cur_site).raw_data(tnd).nbr_frames;
%         if length(nbr_frames)>5
%             for xnd=1:length(nbr_frames)
%                 if nbr_frames(xnd)==5000
%                     if xnd==length(nbr_frames)
%                         order{tnd}(xnd)=3;
%                     else
%                         order{tnd}(xnd)=1;
%                     end
%                 else
%                     order{tnd}(xnd)=2;
%                 end
%             end
%             ord_ind=1:length(order{tnd});
%             order{tnd}(max(ord_ind(order{tnd}==2))+1:end)=3;
%         else
%             order{tnd}=[1 2 1 2 3];
%         end
%         exp_inds{tnd}(1,1)=1;
%         exp_inds{tnd}(1,2)=nbr_frames(1);
%         for fnd=2:length(order{tnd})
%             exp_inds{tnd}(fnd,1)=exp_inds{tnd}(fnd-1,2)+1;
%             exp_inds{tnd}(fnd,2)=sum(nbr_frames(1:fnd));
%         end
%     end
%     vis_fl_all={};
%     delta_all={};
%     for tnd=1:7
%         stim_id=proj_meta(cur_site).raw_data(tnd).stim_id;
%         runM=proj_meta(cur_site).raw_data(tnd).runM;
%         runP=proj_meta(cur_site).raw_data(tnd).runP;
%         pupil_pos=proj_meta(cur_site).raw_data(tnd).pupil_pos;
%         vis_fl=[];
%         cnt=1;
%         for bnd=1:length(order{tnd})
%             if order{tnd}(bnd)==1
%                 tmp_runM=runM(exp_inds{tnd}(bnd,1):exp_inds{tnd}(bnd,2));
%                 tmp_stim_id=stim_id(exp_inds{tnd}(bnd,1):exp_inds{tnd}(bnd,2));
%                 for fnd=1:length(tmp_runM)
%                     if tmp_runM(fnd)==0
%                         vis_fl(cnt)=0;
%                         cnt=cnt+1;
%                     elseif tmp_runM(fnd)==1 & logical(tmp_stim_id(fnd))==0 %feedback
%                         vis_fl(cnt)=1;
%                         cnt=cnt+1;
%                     elseif tmp_runM(fnd)==1 & logical(tmp_stim_id(fnd))==1 %feedback mismatch
%                         if tmp_stim_id(fnd)==7
%                             vis_fl(cnt)=0;
%                             cnt=cnt+1;
%                         else
%                             vis_fl(cnt)=1;
%                             cnt=cnt+1;
%                         end
%                     end
%                 end
%             elseif order{tnd}(bnd)==2
%                 tmp_runM=runM(exp_inds{tnd}(bnd,1):exp_inds{tnd}(bnd,2));
%                 tmp_stim_id=stim_id(exp_inds{tnd}(bnd,1):exp_inds{tnd}(bnd,2));
%                 tmp_runP=runP(exp_inds{tnd}(bnd,1):exp_inds{tnd}(bnd,2));
%                 for fnd=1:length(tmp_runM)
%                     if logical(tmp_stim_id(fnd))==0 & tmp_runM(fnd)==0 & tmp_runP(fnd)==0 | logical(tmp_stim_id(fnd))==1 & tmp_runM(fnd)==0 & tmp_runP(fnd)==0
%                         vis_fl(cnt)=0;
%                         cnt=cnt+1;
%                     elseif logical(tmp_stim_id(fnd))==1 & tmp_runM(fnd)==0 & tmp_runP(fnd)==1
%                         if tmp_stim_id(fnd)==7
%                             vis_fl(cnt)=0;
%                             cnt=cnt+1;
%                         else
%                             vis_fl(cnt)=1;
%                             cnt=cnt+1;
%                         end
%                     elseif logical(tmp_stim_id(fnd))==0 & tmp_runM(fnd)==0 & tmp_runP(fnd)==1
%                         vis_fl(cnt)=1;
%                         cnt=cnt+1;
%                     elseif logical(tmp_stim_id(fnd))==0 & tmp_runM(fnd)==1 & tmp_runP(fnd)==1
%                         vis_fl(cnt)=1;
%                         cnt=cnt+1;
%                     elseif logical(tmp_stim_id(fnd))==0 & tmp_runM(fnd)==1 & tmp_runP(fnd)==0 | logical(tmp_stim_id(fnd))==1 & tmp_runM(fnd)==1 & tmp_runP(fnd)==0
%                         vis_fl(cnt)=0;
%                         cnt=cnt+1;
%                     elseif logical(tmp_stim_id(fnd))==1 & tmp_runM(fnd)==1 & tmp_runP(fnd)==1
%                         if tmp_stim_id(fnd)==7
%                             vis_fl(cnt)=0;
%                             cnt=cnt+1;
%                         else
%                             vis_fl(cnt)=1;
%                             cnt=cnt+1;
%                         end
%                     end
%                 end
%             end
%         end
%         vis_fl(length(vis_fl):length(runM))=0;
%         vis_fl_all{tnd}=vis_fl;
%         delta=[];
%         for ind=1:size(pupil_pos,2)-1
%             dx=abs(pupil_pos(1,ind)-pupil_pos(1,ind+1));
%             dy=abs(pupil_pos(2,ind)-pupil_pos(2,ind+1));
%             delta(ind)=sqrt(dx^2+dy^2);
%             if ind==size(pupil_pos,2)-1
%                 delta(ind+1)=delta(ind);
%             end
%         end
%         delta_all{tnd}=delta;
%     end
%     for fnd=1:size(proj_meta(cur_site).raw_data(1).sig_act,2)
%         %         cur_cell=fnd;
%         snips={};
%         for tnd=1:7
%             raw_act=piezo_act(lnd).activity{tnd}(fnd,:);
%             %             raw_act=all_act{tnd}(fnd,:);
%             cnt=1;
%             ons=find(diff(vis_fl_all{tnd})==1);
%             for gnd=1:length(ons)
%                 if ons(gnd)-low_t>0 && ons(gnd)+high_t<length(vis_fl_all{tnd})
%                     cur_low=ons(gnd)-low_t:ons(gnd);
%                     cur_high=ons(gnd)+1:ons(gnd)+high_t;
%                     if sum(vis_fl_all{tnd}(cur_low))==0 && sum(vis_fl_all{tnd}(cur_high))==high_t
%                         snips{tnd}(cnt,:)=raw_act([cur_low cur_high]);
%                         cnt=cnt+1;
%                     end
%                 end
%             end
%         end
%         exp_av=[];
%         for tnd=1:7
%             if ~isempty(snips{tnd})
%                 exp_av(tnd,:)=mean(snips{tnd},1);
%             else
%                 exp_av(tnd,:)=0;
%             end
%         end
%         if sum(mean(exp_av(:,1:3),2)*1.1<mean(exp_av(:,snip_length-10:snip_length-8),2))
%             res(cnt2,1)=lnd;
%             res(cnt2,2)=fnd;
%             cnt2=cnt2+1;
%         end
%     end
% end
% 
% %take out any runM cells out of vis cells
% vis_res=res;
% remove=[];
% for snd=1:size(run_res,1)
%     if sum(sum(ismember(vis_res,run_res(snd,:)),2)==2)
%         remove(snd)=snd;
%     end
% end
% 
% 
% low_t=5;
% high_t=20;
% snip_length=low_t+high_t+1;
% all_res={};
% for nnd=1:size(res,1)
%     snips={};
%     for tnd=1:7
%         raw_act=piezo_act(res(nnd,1)).activity{tnd}(res(nnd,2),:);
%         cnt=1;
%         ons=find(diff(vis_fl_all{tnd})==1);
% %         vis_ons=vis_ons;
%         for gnd=1:length(ons)
%             if ons(gnd)-low_t>0 && ons(gnd)+high_t<length(vis_fl_all{tnd})
%                 cur_low=ons(gnd)-low_t:ons(gnd);
%                 cur_high=ons(gnd)+1:ons(gnd)+high_t;
%                 if sum(vis_fl_all{tnd}(cur_low))==0 && sum(vis_fl_all{tnd}(cur_high))==high_t
%                     snips{tnd}(cnt,:)=raw_act([cur_low cur_high]);
%                     cnt=cnt+1;
%                 end
%             end
%         end
%     end
%     exp_av=[];
%     for tnd=1:7
% %         if ~isempty(snips{tnd})
% %             exp_av(tnd,:)=
% %         else
% %             exp_av(tnd,:)=0;
% %         end
%         all_res{tnd}(nnd,:)=mean(snips{tnd},1);
%     end
% end
% 
% vis_all_res=all_res;
% figure;
% for tnd=1:7
%     set(gcf,'position',[2003         649        1800         336])
%     subplot(1,7,tnd)
%     imagesc(vis_all_res{tnd})
%     set(gca,'clim',[0.8 2.1])
% end
% title('visual flow onset')
% 
